{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "execution": {},
    "id": "view-in-github"
   },
   "source": [
    "<a href=\"https://colab.research.google.com/github/NeuromatchAcademy/course-content/blob/master/tutorials/W1D3_ModelFitting/student/W1D3_Tutorial1.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a> &nbsp; <a href=\"https://kaggle.com/kernels/welcome?src=https://raw.githubusercontent.com/NeoNeuron/professional-workshop-3/master/tutorials/W4_ModelFitting/student/W4_Tutorial1.ipynb\" target=\"_parent\"><img src=\"https://kaggle.com/static/images/open-in-kaggle.svg\" alt=\"Open in Kaggle\"/></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "execution": {}
   },
   "source": [
    "# Model Fitting: Linear regression with MSE\n",
    "\n",
    "**Original content creators**: Pierre-Étienne Fiquet, Anqi Wu, Alex Hyafil with help from Byron Galbraith\n",
    "\n",
    "**Content modified**: Kai Chen\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "execution": {}
   },
   "source": [
    "___\n",
    "# Tutorial Objectives\n",
    "\n",
    "We start with simple linear regression, using least squares optimization, and then extend to multiple linear regression (Tutorial 1).\n",
    "In the bonus tutorials, we will an alternative methods, Maximum Likelihood Estimation, for linear regression problem, and try to compare the difference with previous ones. (Tutorial Bonus 1).\n",
    "We will use bootstrapping to build confidence intervals around the inferred linear model parameters (Tutorial Bonus 2).\n",
    "We end by learning how to choose between these various models.\n",
    "We discuss the bias-variance trade-off (Tutorial Bonus 3) and Cross Validation for model selection (Tutorial Bonus 4).\n",
    "\n",
    "In this tutorial, we will learn how to fit linear models to data.\n",
    "- Learn how to calculate the mean-squared error (MSE) \n",
    "- Explore how model parameters (slope) influence the MSE\n",
    "- Learn how to find the optimal model parameter using least-squares optimization\n",
    "- Learn how to structure inputs for regression using the 'Design Matrix'\n",
    "- Generalize the MSE for multiple features using the ordinary least squares estimator\n",
    "- Visualize data and model fit in multiple dimensions\n",
    "\n",
    "- Fit polynomial regression models of different complexity\n",
    "- Plot and evaluate the polynomial regression fits\n",
    "\n",
    "---\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "execution": {}
   },
   "source": [
    "---\n",
    "# Setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "both",
    "colab_type": "code",
    "execution": {}
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab_type": "code",
    "execution": {}
   },
   "outputs": [],
   "source": [
    "#@title Figure Settings\n",
    "import ipywidgets as widgets       # interactive display\n",
    "%config InlineBackend.figure_format = 'retina'\n",
    "\n",
    "nma_style = {\n",
    "    'figure.figsize' : (8, 6),\n",
    "    'figure.autolayout' : True,\n",
    "    'font.size' : 15,\n",
    "    'xtick.labelsize' : 'small',\n",
    "    'ytick.labelsize' : 'small',\n",
    "    'legend.fontsize' : 'small',\n",
    "    'axes.spines.top' : False,\n",
    "    'axes.spines.right' : False,\n",
    "    'xtick.major.size' : 5,\n",
    "    'ytick.major.size' : 5,\n",
    "}\n",
    "for key, value in nma_style.items():\n",
    "    plt.rcParams[key] = value\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab_type": "code",
    "execution": {}
   },
   "outputs": [],
   "source": [
    "#@title Helper functions\n",
    "\n",
    "def plot_observed_vs_predicted(x, y, y_hat, theta_hat):\n",
    "  \"\"\" Plot observed vs predicted data\n",
    "\n",
    "  Args:\n",
    "      x (ndarray): observed x values\n",
    "  y (ndarray): observed y values\n",
    "  y_hat (ndarray): predicted y values\n",
    "\n",
    "  \"\"\"\n",
    "  fig, ax = plt.subplots()\n",
    "  ax.scatter(x, y, label='Observed')  # our data scatter plot\n",
    "  ax.plot(x, y_hat, color='r', label='Fit')  # our estimated model\n",
    "  # plot residuals\n",
    "  ymin = np.minimum(y, y_hat)\n",
    "  ymax = np.maximum(y, y_hat)\n",
    "  ax.vlines(x, ymin, ymax, 'g', alpha=0.5, label='Residuals')\n",
    "  ax.set(\n",
    "      title=fr\"$\\hat{{\\theta}}$ = {theta_hat:0.2f}, MSE = {mse(x, y, theta_hat):.2f}\",\n",
    "      xlabel='x',\n",
    "      ylabel='y'\n",
    "  )\n",
    "  ax.legend()\n",
    "\n",
    "def plot_fitted_polynomials(x, y, theta_hat):\n",
    "  \"\"\" Plot polynomials of different orders\n",
    "\n",
    "  Args:\n",
    "    x (ndarray): input vector of shape (n_samples)\n",
    "    y (ndarray): vector of measurements of shape (n_samples)\n",
    "    theta_hat (dict): polynomial regression weights for different orders\n",
    "  \"\"\"\n",
    "\n",
    "  x_grid = np.linspace(x.min() - .5, x.max() + .5)\n",
    "\n",
    "  plt.figure()\n",
    "\n",
    "  for order in range(0, max_order + 1):\n",
    "    X_design = make_design_matrix(x_grid, order)\n",
    "    plt.plot(x_grid, X_design @ theta_hat[order]);\n",
    "\n",
    "  plt.ylabel('y')\n",
    "  plt.xlabel('x')\n",
    "  plt.plot(x, y, 'C0.');\n",
    "  plt.legend([f'order {o}' for o in range(max_order + 1)], loc=1)\n",
    "  plt.title('polynomial fits')\n",
    "  plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "execution": {}
   },
   "source": [
    "---\n",
    "# Section 1: Mean Squared Error (MSE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "execution": {}
   },
   "source": [
    "**Linear least squares regression** is an old but gold  optimization procedure that we are going to use for data fitting. Least squares (LS) optimization problems are those in which the objective function is a quadratic function of the\n",
    "parameter(s) being optimized.\n",
    "\n",
    "Suppose you have a set of measurements, $y_{n}$ (the \"dependent\" variable) obtained for different input values, $x_{n}$ (the \"independent\" or \"explanatory\" variable). Suppose we believe the measurements are proportional to the input values, but are corrupted by some (random) measurement errors, $\\epsilon_{n}$, that is:\n",
    "\n",
    "$$y_{n}= \\theta x_{n}+\\epsilon_{n}$$\n",
    "\n",
    "for some unknown slope parameter $\\theta.$ The least squares regression problem uses **mean squared error (MSE)** as its objective function, it aims to find the value of the parameter $\\theta$ by minimizing the average of squared errors:\n",
    "\n",
    "\\begin{align}\n",
    "\\min _{\\theta} \\frac{1}{N}\\sum_{n=1}^{N}\\left(y_{n}-\\theta x_{n}\\right)^{2}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "execution": {}
   },
   "source": [
    "We will now explore how MSE is used in fitting a linear regression model to data. For illustrative purposes, we will create a simple synthetic dataset where we know the true underlying model. This will allow us to see how our estimation efforts compare in uncovering the real model (though in practice we rarely have this luxury).\n",
    "\n",
    "First we will generate some noisy samples $x$ from [0, 10) along the line $y = 1.2x$ as our dataset we wish to fit a model to."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab_type": "code",
    "execution": {}
   },
   "outputs": [],
   "source": [
    "# @title\n",
    "\n",
    "# @markdown Execute this cell to generate some simulated data\n",
    "\n",
    "# setting a fixed seed to our random number generator ensures we will always\n",
    "# get the same psuedorandom number sequence\n",
    "np.random.seed(121)\n",
    "\n",
    "# Let's set some parameters\n",
    "theta = 1.2\n",
    "n_samples = 30\n",
    "\n",
    "# Draw x and then calculate y\n",
    "x = 10 * np.random.rand(n_samples)  # sample from a uniform distribution over [0,10)\n",
    "noise = np.random.randn(n_samples)  # sample from a standard normal distribution\n",
    "y = theta * x + noise\n",
    "\n",
    "# Plot the results\n",
    "fig, ax = plt.subplots()\n",
    "ax.scatter(x, y)  # produces a scatter plot\n",
    "ax.set(xlabel='x', ylabel='y');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "execution": {}
   },
   "source": [
    "Now that we have our suitably noisy dataset, we can start trying to estimate the underlying model that produced it. We use MSE to evaluate how successful a particular slope estimate $\\hat{\\theta}$ is for explaining the data, with the closer to 0 the MSE is, the better our estimate fits the data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "execution": {}
   },
   "source": [
    "## Exercise 1: Compute MSE\n",
    "\n",
    "In this exercise you will implement a method to compute the mean squared error for a set of inputs $x$, measurements $y$, and slope estimate $\\hat{\\theta}$. We will then compute and print the mean squared error for 3 different choices of theta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab_type": "code",
    "execution": {}
   },
   "outputs": [],
   "source": [
    "def mse(x, y, theta_hat):\n",
    "  \"\"\"Compute the mean squared error\n",
    "\n",
    "  Args:\n",
    "    x (ndarray): An array of shape (samples,) that contains the input values.\n",
    "    y (ndarray): An array of shape (samples,) that contains the corresponding\n",
    "      measurement values to the inputs.\n",
    "    theta_hat (float): An estimate of the slope parameter\n",
    "\n",
    "  Returns:\n",
    "    float: The mean squared error of the data with the estimated parameter.\n",
    "  \"\"\"\n",
    "  ####################################################\n",
    "  ## TODO for students: compute the mean squared error\n",
    "  # Fill out function and remove\n",
    "  raise NotImplementedError(\"Student exercise: compute the mean squared error\")\n",
    "  ####################################################\n",
    "\n",
    "  # Compute the estimated y\n",
    "  y_hat = ...\n",
    "\n",
    "  # Compute mean squared error\n",
    "  mse = ...\n",
    "\n",
    "  return mse\n",
    "\n",
    "\n",
    "# Uncomment below to test your function\n",
    "theta_hats = [0.75, 1.0, 1.5]\n",
    "# for theta_hat in theta_hats:\n",
    "#   print(f\"theta_hat of {theta_hat} has an MSE of {mse(x, y, theta_hat):.2f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "cellView": "both",
    "colab_type": "text",
    "execution": {}
   },
   "source": [
    "[*Click for solution*](https://github.com/NeoNeuron/professional-workshop-3/tree/master//tutorials/W4_ModelFitting/solutions/W4_Tutorial1_Solution_12a57de0.py)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "execution": {}
   },
   "source": [
    "The result should be:\n",
    "\n",
    "theta_hat of 0.75 has an MSE of 9.08\\\n",
    "theta_hat of 1.0 has an MSE of 3.0\\\n",
    "theta_hat of 1.5 has an MSE of 4.52\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "execution": {}
   },
   "source": [
    "We see that $\\hat{\\theta} = 1.0$ is our best estimate from the three we tried. Looking just at the raw numbers, however, isn't always satisfying, so let's visualize what our estimated model looks like over the data. \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab_type": "code",
    "execution": {}
   },
   "outputs": [],
   "source": [
    "#@title\n",
    "\n",
    "#@markdown Execute this cell to visualize estimated models\n",
    "\n",
    "fig, axes = plt.subplots(ncols=3, figsize=(18, 4))\n",
    "for theta_hat, ax in zip(theta_hats, axes):\n",
    "\n",
    "  # True data\n",
    "  ax.scatter(x, y, label='Observed')  # our data scatter plot\n",
    "\n",
    "  # Compute and plot predictions\n",
    "  y_hat = theta_hat * x\n",
    "  ax.plot(x, y_hat, color='r', label='Fit')  # our estimated model\n",
    "\n",
    "  ax.set(\n",
    "      title= fr'$\\hat{{\\theta}}$= {theta_hat}, MSE = {mse(x, y, theta_hat):.2f}',\n",
    "      xlabel='x',\n",
    "      ylabel='y'\n",
    "  );\n",
    "\n",
    "axes[0].legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "execution": {}
   },
   "source": [
    "## Interactive Demo: MSE Explorer\n",
    "\n",
    "Using an interactive widget, we can easily see how changing our slope estimate changes our model fit. We display the **residuals**, the differences between observed and predicted data, as line segments between the data point (observed response) and the corresponding predicted response on the model fit line."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab_type": "code",
    "execution": {}
   },
   "outputs": [],
   "source": [
    "#@title\n",
    "\n",
    "#@markdown Make sure you execute this cell to enable the widget!\n",
    "\n",
    "@widgets.interact(theta_hat=widgets.FloatSlider(1.0, min=0.0, max=2.0))\n",
    "def plot_data_estimate(theta_hat):\n",
    "  y_hat = theta_hat * x\n",
    "  plot_observed_vs_predicted(x, y, y_hat, theta_hat)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "execution": {}
   },
   "source": [
    "While visually exploring several estimates can be instructive, it's not the most efficient for finding the best estimate to fit our data. Another technique we can use is choose a reasonable range of parameter values and compute the MSE at several values in that interval. This allows us to plot the error against the parameter value (this is also called an **error landscape**, especially when we deal with more than one parameter). We can select the final $\\hat{\\theta}$ ($\\hat{\\theta}_{MSE}$) as the one which results in the lowest error."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab_type": "code",
    "execution": {}
   },
   "outputs": [],
   "source": [
    "# @title\n",
    "\n",
    "# @markdown Execute this cell to loop over theta_hats, compute MSE, and plot results\n",
    "\n",
    "# Loop over different thetas, compute MSE for each\n",
    "theta_hat_grid = np.linspace(-2.0, 4.0)\n",
    "errors = np.zeros(len(theta_hat_grid))\n",
    "for i, theta_hat in enumerate(theta_hat_grid):\n",
    "  errors[i] = mse(x, y, theta_hat)\n",
    "\n",
    "# Find theta that results in lowest error\n",
    "best_error = np.min(errors)\n",
    "theta_hat = theta_hat_grid[np.argmin(errors)]\n",
    "\n",
    "\n",
    "# Plot results\n",
    "fig, ax = plt.subplots()\n",
    "ax.plot(theta_hat_grid, errors, '-o', label='MSE', c='C1')\n",
    "ax.axvline(theta, color='g', ls='--', label=r\"$\\theta_{True}$\")\n",
    "ax.axvline(theta_hat, color='r', ls='-', label=r\"$\\hat{{\\theta}}_{MSE}$\")\n",
    "ax.set(\n",
    "  title=fr\"Best fit: $\\hat{{\\theta}}$ = {theta_hat:.2f}, MSE = {best_error:.2f}\",\n",
    "  xlabel=r\"$\\hat{{\\theta}}$\",\n",
    "  ylabel='MSE')\n",
    "ax.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "execution": {}
   },
   "source": [
    "We can see that our best fit is $\\hat{\\theta}=1.18$ with an MSE of 1.45. This is quite close to the original true value $\\theta=1.2$!\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "execution": {}
   },
   "source": [
    "---\n",
    "# Section 2: Least-squares optimization\n",
    "\n",
    "\n",
    "While the approach detailed above (computing MSE at various values of $\\hat\\theta$) quickly got us to a good estimate, it still relied on evaluating the MSE value across a grid of hand-specified values. If we didn't pick a good range to begin with, or with enough granularity, we might miss the best possible estimator. Let's go one step further, and instead of finding the minimum MSE from a set of candidate estimates, let's solve for it analytically.\n",
    "\n",
    "We can do this by minimizing the cost function. Mean squared error is a convex objective function, therefore we can compute its minimum using calculus. Please see video or appendix for this derivation! After computing the minimum, we find that:\n",
    "\n",
    "\\begin{align}\n",
    "\\hat\\theta = \\frac{\\vec{x}^\\top \\vec{y}}{\\vec{x}^\\top \\vec{x}}\n",
    "\\end{align}\n",
    "\n",
    "This is known as solving the normal equations. For different ways of obtaining the solution, see the notes on [Least Squares Optimization](https://www.cns.nyu.edu/~eero/NOTES/leastSquares.pdf) by Eero Simoncelli."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "execution": {}
   },
   "source": [
    "### Exercise 2: Solve for the Optimal Estimator\n",
    "\n",
    "In this exercise, you will write a function that finds the optimal $\\hat{\\theta}$ value using the least squares optimization approach (the equation above) to solve MSE minimization. It shoud take arguments $x$ and $y$ and return the solution $\\hat{\\theta}$.\n",
    "\n",
    "We will then use your function to compute $\\hat{\\theta}$ and plot the resulting prediction on top of the data. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab_type": "code",
    "execution": {}
   },
   "outputs": [],
   "source": [
    "def solve_normal_eqn(x, y):\n",
    "  \"\"\"Solve the normal equations to produce the value of theta_hat that minimizes\n",
    "    MSE.\n",
    "\n",
    "    Args:\n",
    "    x (ndarray): An array of shape (samples,) that contains the input values.\n",
    "    y (ndarray): An array of shape (samples,) that contains the corresponding\n",
    "      measurement values to the inputs.\n",
    "\n",
    "  Returns:\n",
    "    float: the value for theta_hat arrived from minimizing MSE\n",
    "  \"\"\"\n",
    "\n",
    "  ################################################################################\n",
    "  ## TODO for students: solve for the best parameter using least squares\n",
    "  # Fill out function and remove\n",
    "  raise NotImplementedError(\"Student exercise: solve for theta_hat using least squares\")\n",
    "  ################################################################################\n",
    "\n",
    "  # Compute theta_hat analytically\n",
    "  theta_hat = ...\n",
    "\n",
    "  return theta_hat\n",
    "\n",
    "\n",
    "# Uncomment below to test your function\n",
    "# theta_hat = solve_normal_eqn(x, y)\n",
    "# y_hat = theta_hat * x\n",
    "# plot_observed_vs_predicted(x, y, y_hat, theta_hat)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "cellView": "both",
    "colab_type": "text",
    "execution": {}
   },
   "source": [
    "[*Click for solution*](https://github.com/NeoNeuron/professional-workshop-3/tree/master//tutorials/W4_ModelFitting/solutions/W4_Tutorial1_Solution_b1fcbf07.py)\n",
    "\n",
    "*Example output:*\n",
    "\n",
    "<img alt='Solution hint' align='left' width=1114.0 height=827.0 src=https://raw.githubusercontent.com/NeoNeuron/professional-workshop-3/master/tutorials/W4_ModelFitting/static/W4_Tutorial1_Solution_b1fcbf07_0.png>\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "execution": {}
   },
   "source": [
    "We see that the analytic solution produces an even better result than our grid search from before, producing $\\hat{\\theta} = 1.21$ with MSE = 1.43!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "---\n",
    "# Section 3: Multiple Linear Regression\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "\n",
    "Now that we have considered the univariate case and how to produce confidence intervals for our estimator, we turn to the general linear regression case, where we can have more than one regressor, or feature, in our input.\n",
    "\n",
    "Recall that our original univariate linear model was given as\n",
    "\n",
    "\\begin{align}\n",
    "y = \\theta x + \\epsilon\n",
    "\\end{align}\n",
    "\n",
    "where $\\theta$ is the slope and $\\epsilon$ some noise. We can easily extend this to the multivariate scenario by adding another parameter for each additional feature\n",
    "\n",
    "\\begin{align}\n",
    "y = \\theta_0 + \\theta_1 x_1 + \\theta_1 x_2 + ... +\\theta_d x_d + \\epsilon\n",
    "\\end{align}\n",
    "\n",
    "where $\\theta_0$ is the intercept and $d$ is the number of features (it is also the dimensionality of our input).\n",
    "\n",
    "We can condense this succinctly using vector notation for a single data point\n",
    "\n",
    "\\begin{align}\n",
    "y_i = \\boldsymbol{\\theta}^{\\top}\\mathbf{x}_i + \\epsilon\n",
    "\\end{align}\n",
    "\n",
    "and fully in matrix form\n",
    "\n",
    "\\begin{align}\n",
    "\\mathbf{y} = \\mathbf{X}\\boldsymbol{\\theta} + \\mathbf{\\epsilon}\n",
    "\\end{align}\n",
    "\n",
    "where $\\mathbf{y}$ is a vector of measurements, $\\mathbf{X}$ is a matrix containing the feature values (columns) for each input sample (rows), and $\\boldsymbol{\\theta}$ is our parameter vector.\n",
    "\n",
    "This matrix $\\mathbf{X}$ is often referred to as the \"[design matrix](https://en.wikipedia.org/wiki/Design_matrix)\"."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "\n",
    "Here we will focus on the two-dimensional case ($d=2$), which allows us to easily visualize our results. As an example, think of a situation where a scientist records the spiking response of a retinal ganglion cell to patterns of light signals that vary in contrast and in orientation. Then contrast and orientation values can be used as features / regressors to predict the cells response.\n",
    "\n",
    "In this case our model can be writen as:\n",
    "\n",
    "\\begin{align}\n",
    "y = \\theta_0 + \\theta_1 x_1 + \\theta_2 x_2 + \\epsilon\n",
    "\\end{align}\n",
    "\n",
    "or in matrix form where\n",
    "\n",
    "\\begin{align}\n",
    "\\mathbf{X} = \n",
    "\\begin{bmatrix}\n",
    "1 & x_{1,1} & x_{1,2} \\\\\n",
    "1 & x_{2,1} & x_{2,2} \\\\\n",
    "\\vdots & \\vdots & \\vdots \\\\\n",
    "1 & x_{n,1} & x_{n,2}\n",
    "\\end{bmatrix}, \n",
    "\\boldsymbol{\\theta} =\n",
    "\\begin{bmatrix}\n",
    "\\theta_0 \\\\\n",
    "\\theta_1 \\\\\n",
    "\\theta_2 \\\\\n",
    "\\end{bmatrix}\n",
    "\\end{align}\n",
    "\n",
    "For our actual exploration dataset we shall set $\\boldsymbol{\\theta}=[0, -2, -3]$ and draw $N=40$ noisy samples from $x \\in [-2,2)$. Note that setting the value of $\\theta_0 = 0$ effectively ignores the offset term."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "execution": {}
   },
   "outputs": [],
   "source": [
    "#@title\n",
    "#@markdown Execute this cell to simulate some data\n",
    "\n",
    "# Set random seed for reproducibility\n",
    "np.random.seed(1234)\n",
    "\n",
    "# Set parameters\n",
    "theta = [0, -2, -3]\n",
    "n_samples = 40\n",
    "\n",
    "# Draw x and calculate y\n",
    "n_regressors = len(theta)\n",
    "x0 = np.ones((n_samples, 1))\n",
    "x1 = np.random.uniform(-2, 2, (n_samples, 1))\n",
    "x2 = np.random.uniform(-2, 2, (n_samples, 1))\n",
    "X = np.hstack((x0, x1, x2))\n",
    "noise = np.random.randn(n_samples)\n",
    "y = X @ theta + noise\n",
    "\n",
    "\n",
    "ax = plt.subplot(projection='3d')\n",
    "ax.plot(X[:,1], X[:,2], y, '.')\n",
    "\n",
    "ax.set(\n",
    "    xlabel='$x_1$',\n",
    "    ylabel='$x_2$',\n",
    "    zlabel='y'\n",
    ")\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "Now that we have our dataset, we want to find an optimal vector of paramters $\\boldsymbol{\\hat\\theta}$. Recall our analytic solution to minimizing MSE for a single regressor:\n",
    "\n",
    "\\begin{align}\n",
    "\\hat\\theta = \\frac{\\sum_{i=1}^N x_i y_i}{\\sum_{i=1}^N x_i^2}.\n",
    "\\end{align}\n",
    "\n",
    "The same holds true for the multiple regressor case, only now expressed in matrix form\n",
    "\n",
    "\\begin{align}\n",
    "\\boldsymbol{\\hat\\theta} = (\\mathbf{X}^\\top\\mathbf{X})^{-1}\\mathbf{X}^\\top\\mathbf{y}.\n",
    "\\end{align}\n",
    "\n",
    "This is called the [ordinary least squares](https://en.wikipedia.org/wiki/Ordinary_least_squares) (OLS) estimator."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "### Exercise 3: Ordinary Least Squares Estimator\n",
    "\n",
    "In this exercise you will implement the OLS approach to estimating $\\boldsymbol{\\hat\\theta}$ from the design matrix $\\mathbf{X}$ and measurement vector $\\mathbf{y}$. You can use the `@` symbol for matrix multiplication, `.T` for transpose, and `np.linalg.inv` for matrix inversion.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {}
   },
   "outputs": [],
   "source": [
    "def ordinary_least_squares(X, y):\n",
    "  \"\"\"Ordinary least squares estimator for linear regression.\n",
    "\n",
    "  Args:\n",
    "    x (ndarray): design matrix of shape (n_samples, n_regressors)\n",
    "    y (ndarray): vector of measurements of shape (n_samples)\n",
    "\n",
    "  Returns:\n",
    "    ndarray: estimated parameter values of shape (n_regressors)\n",
    "  \"\"\"\n",
    "  ######################################################################\n",
    "  ## TODO for students: solve for the optimal parameter vector using OLS\n",
    "  # Fill out function and remove\n",
    "  raise NotImplementedError(\"Student exercise: solve for theta_hat vector using OLS\")\n",
    "  ######################################################################\n",
    "\n",
    "  # Compute theta_hat using OLS\n",
    "  theta_hat = ...\n",
    "\n",
    "  return theta_hat\n",
    "\n",
    "\n",
    "# Uncomment below to test your function\n",
    "# theta_hat = ordinary_least_squares(X, y)\n",
    "# print(theta_hat)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "execution": {}
   },
   "source": [
    "[*Click for solution*](https://github.com/NeoNeuron/professional-workshop-3/tree/master//tutorials/W4_ModelFitting/solutions/W4_Tutorial1_Solution_2b5abc38.py)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "After filling in this function, you should see that $\\hat{\\theta}$ = [ 0.13861386, -2.09395731, -3.16370742]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "Now that we have our $\\mathbf{\\hat\\theta}$, we can obtain $\\mathbf{\\hat y}$ and thus our mean squared error."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {}
   },
   "outputs": [],
   "source": [
    "# Compute predicted data\n",
    "theta_hat = ordinary_least_squares(X, y)\n",
    "y_hat = X @ theta_hat\n",
    "\n",
    "# Compute MSE\n",
    "print(f\"MSE = {np.mean((y - y_hat)**2):.2f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "Finally, the following code will plot a geometric visualization of the data points (blue) and fitted plane. You can see that the residuals (green bars) are orthogonal to the plane.\n",
    "\n",
    "Extra: how would you check that the residuals are orthogonal to the fitted plane? (this is sometimes called the orthogonality principle, see [least squares notes](https://www.cns.nyu.edu/~eero/NOTES/leastSquares.pdf))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "execution": {}
   },
   "outputs": [],
   "source": [
    "#@title\n",
    "#@markdown Execute this cell to visualize data and predicted plane\n",
    "\n",
    "theta_hat = ordinary_least_squares(X, y)\n",
    "xx, yy = np.mgrid[-2:2:50j, -2:2:50j]\n",
    "y_hat_grid = np.array([xx.flatten(), yy.flatten()]).T @ theta_hat[1:]\n",
    "y_hat_grid = y_hat_grid.reshape((50, 50))\n",
    "\n",
    "ax = plt.subplot(projection='3d')\n",
    "ax.plot(X[:, 1], X[:, 2], y, '.')\n",
    "ax.plot_surface(xx, yy, y_hat_grid, linewidth=0, alpha=0.5, color='C1',\n",
    "                cmap=plt.get_cmap('coolwarm'))\n",
    "\n",
    "for i in range(len(X)):\n",
    "  ax.plot((X[i, 1], X[i, 1]),\n",
    "          (X[i, 2], X[i, 2]),\n",
    "          (y[i], y_hat[i]),\n",
    "          'g-', alpha=.5)\n",
    "\n",
    "ax.set(\n",
    "    xlabel='$x_1$',\n",
    "    ylabel='$x_2$',\n",
    "    zlabel='y'\n",
    ")\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "execution": {}
   },
   "source": [
    "---\n",
    "# Summary\n",
    "\n",
    "- Linear least squares regression is an optimization procedure that can be used for data fitting:\n",
    "    - Task: predict a value for $y$ given $x$\n",
    "    - Performance measure: $\\textrm{MSE}$\n",
    "    - Procedure: minimize $\\textrm{MSE}$ by solving the normal equations\n",
    "- Linear regression generalizes naturally to multiple dimensions\n",
    "- Linear algebra affords us the mathematical tools to reason and solve such problems beyond the two dimensional case\n",
    "- **Key point**: We fit the model by defining an *objective function* and minimizing it. \n",
    "- **Note**: \n",
    "    - In this case, there is an *analytical* solution to the minimization problem and in practice, this solution can be computed using *linear algebra*. This is *extremely* powerful and forms the basis for much of numerical computation throughout the sciences.\n",
    "    - In practice, multidimensional least squares problems can be solved very efficiently (thanks to numerical routines such as LAPACK)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "execution": {}
   },
   "source": [
    "---\n",
    "# Bonus:  Polynomial Regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "execution": {}
   },
   "source": [
    "So far today, you learned how to predict outputs from inputs by fitting a linear regression model. We can now model all sort of relationships, including in neuroscience!  \n",
    "\n",
    "One potential problem with this approach is the simplicity of the model. Linear regression, as the name implies, can only capture a linear relationship between the inputs and outputs. Put another way, the predicted outputs are only a weighted sum of the inputs. What if there are more complicated computations happening? Luckily, many more complex models exist (and you will encounter many more over the next 3 weeks). One model that is still very simple to fit and understand, but captures more complex relationships, is **polynomial regression**, an extension of linear regression.\n",
    "\n",
    "Since polynomial regression is an extension of linear regression, everything you learned so far will come in handy now! The goal is the same: we want to predict the dependent variable $y_{n}$ given the input values $x_{n}$. The key change is the type of relationship between inputs and outputs that the model can capture. \n",
    "\n",
    "Linear regression models predict the outputs as a weighted sum of the inputs:\n",
    "\n",
    "$$y_{n}= \\theta_0 + \\theta x_{n} + \\epsilon_{n}$$\n",
    "\n",
    "With polynomial regression, we model the outputs as a polynomial equation based on the inputs. For example, we can model the outputs as:\n",
    "\n",
    "$$y_{n}= \\theta_0 + \\theta_1 x_{n} + \\theta_2 x_{n}^2 + \\theta_3 x_{n}^3 + \\epsilon_{n}$$\n",
    "\n",
    "We can change how complex a polynomial is fit by changing the order of the polynomial. The order of a polynomial refers to the highest power in the polynomial. The equation above is a third order polynomial because the highest value x is raised to is 3. We could add another term ($+ \\theta_4 x_{n}^4$) to model an order 4 polynomial and so on.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "execution": {}
   },
   "source": [
    "First, we will simulate some data to practice fitting polynomial regression models. We will generate random inputs $x$ and then compute y according to $y = x^2 - x - 2 $, with some extra noise both in the input and the output to make the model fitting exercise closer to a real life situation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab_type": "code",
    "execution": {}
   },
   "outputs": [],
   "source": [
    "#@title\n",
    "#@markdown Execute this cell to simulate some data\n",
    "\n",
    "# setting a fixed seed to our random number generator ensures we will always\n",
    "# get the same psuedorandom number sequence\n",
    "np.random.seed(121)\n",
    "n_samples = 30\n",
    "x = np.random.uniform(-2, 2.5, n_samples)  # inputs uniformly sampled from [-2, 2.5)\n",
    "y =  x**2 - x - 2   # computing the outputs\n",
    "\n",
    "output_noise = 1/8 * np.random.randn(n_samples)\n",
    "y += output_noise  # adding some output noise\n",
    "\n",
    "input_noise = 1/2 * np.random.randn(n_samples)\n",
    "x += input_noise  # adding some input noise\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "ax.scatter(x, y)  # produces a scatter plot\n",
    "ax.set(xlabel='x', ylabel='y');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "execution": {}
   },
   "source": [
    "## Design matrix for polynomial regression\n",
    "\n",
    "Now we have the basic idea of polynomial regression and some noisy data, let's begin! The key difference between fitting a linear regression model and a polynomial regression model lies in how we structure the input variables.  \n",
    "\n",
    "For linear regression, we used $X = x$ as the input data. To add a constant bias (a y-intercept in a 2-D plot), we use $X = \\big[ \\boldsymbol 1, x \\big]$, where $\\boldsymbol 1$ is a column of ones.  When fitting, we learn a weight for each column of this matrix. So we learn a weight that multiples with column 1 - in this case that column is all ones so we gain the bias parameter ($+ \\theta_0$). We also learn a weight for every column, or every feature of x, as we learned in Section 1.\n",
    "\n",
    "This matrix $X$ that we use for our inputs is known as a **design matrix**. We want to create our design matrix so we learn weights for $x^2, x^3,$ etc. Thus, we want to build our design matrix $X$ for polynomial regression of order $k$ as:\n",
    "\n",
    "$$X = \\big[ \\boldsymbol 1 , x^1, x^2 , \\ldots , x^k \\big],$$\n",
    "\n",
    "where $\\boldsymbol{1}$ is the vector the same length as $x$ consisting of of all ones, and $x^p$ is the vector or matrix $x$ with all elements raised to the power $p$. Note that $\\boldsymbol{1} = x^0$ and $x^1 = x$ "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "execution": {}
   },
   "source": [
    "### Exercise : Structure design matrix\n",
    "\n",
    "Create a function (`make_design_matrix`) that structures the design matrix given the input data and the order of the polynomial you wish to fit. We will print part of this design matrix for our data and order 5."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab_type": "code",
    "execution": {}
   },
   "outputs": [],
   "source": [
    "def make_design_matrix(x, order):\n",
    "  \"\"\"Create the design matrix of inputs for use in polynomial regression\n",
    "\n",
    "  Args:\n",
    "    x (ndarray): input vector of shape (n_samples)\n",
    "    order (scalar): polynomial regression order\n",
    "\n",
    "  Returns:\n",
    "    ndarray: design matrix for polynomial regression of shape (samples, order+1)\n",
    "  \"\"\"\n",
    "  ########################################################################\n",
    "  ## TODO for students: create the design matrix ##\n",
    "  # Fill out function and remove\n",
    "  raise NotImplementedError(\"Student exercise: create the design matrix\")\n",
    "  ########################################################################\n",
    "\n",
    "  # Broadcast to shape (n x 1) so dimensions work\n",
    "  if x.ndim == 1:\n",
    "    x = x[:, None]\n",
    "\n",
    "  #if x has more than one feature, we don't want multiple columns of ones so we assign\n",
    "  # x^0 here\n",
    "  design_matrix = np.ones((x.shape[0], 1))\n",
    "\n",
    "  # Loop through rest of degrees and stack columns (hint: np.hstack)\n",
    "  for degree in range(1, order + 1):\n",
    "      design_matrix = ...\n",
    "\n",
    "  return design_matrix\n",
    "\n",
    "\n",
    "# Uncomment to test your function\n",
    "order = 5\n",
    "# X_design = make_design_matrix(x, order)\n",
    "# print(X_design[0:2, 0:2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "cellView": "both",
    "colab_type": "text",
    "execution": {}
   },
   "source": [
    "[*Click for solution*](https://github.com/NeoNeuron/professional-workshop-3/tree/master//tutorials/W4_ModelFitting/solutions/W4_Tutorial1_Solution_f8576b48.py)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "execution": {}
   },
   "source": [
    "You should see that the printed section of this design matrix is `[[ 1.         -1.51194917]\n",
    " [ 1.         -0.35259945]]`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "execution": {}
   },
   "source": [
    "## Fitting polynomial regression models"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "execution": {}
   },
   "source": [
    "Now that we have the inputs structured correctly in our design matrix, fitting a polynomial regression is the same as fitting a linear regression model! All of the polynomial structure we need to learn is contained in how the inputs are structured in the design matrix. We can use the same least squares solution we computed in previous exercises. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "execution": {}
   },
   "source": [
    "### Exercise : Fitting polynomial regression models with different orders "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "execution": {}
   },
   "source": [
    "Here, we will fit polynomial regression models to find the regression coefficients ($\\theta_0, \\theta_1, \\theta_2,$ ...) by solving the least squares problem. Create a function `solve_poly_reg` that loops over different order polynomials (up to `max_order`), fits that model, and saves out the weights for each. You may invoke the `ordinary_least_squares` function. \n",
    "\n",
    "We will then qualitatively inspect the quality of our fits for each order by plotting the fitted polynomials on top of the data. In order to see smooth curves, we evaluate the fitted polynomials on a grid of $x$ values (ranging between the largest and smallest of the inputs present in the dataset)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab_type": "code",
    "execution": {}
   },
   "outputs": [],
   "source": [
    "def solve_poly_reg(x, y, max_order):\n",
    "  \"\"\"Fit a polynomial regression model for each order 0 through max_order.\n",
    "\n",
    "  Args:\n",
    "    x (ndarray): input vector of shape (n_samples)\n",
    "    y (ndarray): vector of measurements of shape (n_samples)\n",
    "    max_order (scalar): max order for polynomial fits\n",
    "\n",
    "  Returns:\n",
    "    dict: fitted weights for each polynomial model (dict key is order)\n",
    "  \"\"\"\n",
    "\n",
    "  # Create a dictionary with polynomial order as keys,\n",
    "  # and np array of theta_hat (weights) as the values\n",
    "  theta_hats = {}\n",
    "\n",
    "  # Loop over polynomial orders from 0 through max_order\n",
    "  for order in range(max_order + 1):\n",
    "\n",
    "    ##################################################################################\n",
    "    ## TODO for students: Create design matrix and fit polynomial model for this order\n",
    "    # Fill out function and remove\n",
    "    raise NotImplementedError(\"Student exercise: fit a polynomial model\")\n",
    "    ##################################################################################\n",
    "\n",
    "    # Create design matrix\n",
    "    X_design = ...\n",
    "\n",
    "    # Fit polynomial model\n",
    "    this_theta = ...\n",
    "\n",
    "    theta_hats[order] = this_theta\n",
    "\n",
    "  return theta_hats\n",
    "\n",
    "\n",
    "# Uncomment to test your function\n",
    "max_order = 5\n",
    "# theta_hats = solve_poly_reg(x, y, max_order)\n",
    "# plot_fitted_polynomials(x, y, theta_hats)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "cellView": "both",
    "colab_type": "text",
    "execution": {}
   },
   "source": [
    "[*Click for solution*](https://github.com/NeoNeuron/professional-workshop-3/tree/master//tutorials/W4_ModelFitting/solutions/W4_Tutorial1_Solution_9b1eebd9.py)\n",
    "\n",
    "*Example output:*\n",
    "\n",
    "<img alt='Solution hint' align='left' width=1120.0 height=832.0 src=https://raw.githubusercontent.com/NeoNeuron/professional-workshop-3/master/tutorials/W4_ModelFitting/static/W4_Tutorial1_Solution_9b1eebd9_3.png>\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "execution": {}
   },
   "source": [
    "## Evaluating fit quality\n",
    "\n",
    "As with linear regression, we can compute mean squared error (MSE) to get a sense of how well the model fits the data. \n",
    "\n",
    "We compute MSE as:\n",
    "\n",
    "$$ MSE = \\frac 1 N ||y - \\hat y||^2 = \\sum_{i=1}^N (y_i - \\hat y_i)^2 $$\n",
    "\n",
    "where the predicted values for each model are given by $ \\hat y = X \\hat \\theta$. \n",
    "\n",
    "*Which model (i.e. which polynomial order) do you think will have the best MSE?*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "execution": {}
   },
   "source": [
    "### Exercise : Compute MSE and compare models\n",
    "\n",
    "\n",
    "We will compare the MSE for different polynomial orders with a bar plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab_type": "code",
    "execution": {}
   },
   "outputs": [],
   "source": [
    "mse_list = []\n",
    "order_list = list(range(max_order + 1))\n",
    "\n",
    "for order in order_list:\n",
    "\n",
    "  X_design = make_design_matrix(x, order)\n",
    "\n",
    "  ############################################################################\n",
    "  # TODO for students: Compute MSE (fill in ... sections)\n",
    "  #############################################################################\n",
    "\n",
    "  # Get prediction for the polynomial regression model of this order\n",
    "  y_hat = ...\n",
    "\n",
    "  # Compute the residuals\n",
    "  residuals = ...\n",
    "\n",
    "  # Compute the MSE\n",
    "  mse = ...\n",
    "\n",
    "  mse_list.append(mse)\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "\n",
    "# Uncomment once above exercise is complete\n",
    "# ax.bar(order_list, mse_list)\n",
    "\n",
    "ax.set(title='Comparing Polynomial Fits', xlabel='Polynomial order', ylabel='MSE')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "cellView": "both",
    "colab_type": "text",
    "execution": {}
   },
   "source": [
    "[*Click for solution*](https://github.com/NeoNeuron/professional-workshop-3/tree/master//tutorials/W4_ModelFitting/solutions/W4_Tutorial1_Solution_ccea5dd7.py)\n",
    "\n",
    "*Example output:*\n",
    "\n",
    "<img alt='Solution hint' align='left' width=1116.0 height=827.0 src=https://raw.githubusercontent.com/NeoNeuron/professional-workshop-3/master/tutorials/W4_ModelFitting/static/W4_Tutorial1_Solution_ccea5dd7_0.png>\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "\n",
    "* To change from a linear regression model to a polynomial regression model, we only have to change how the input data is structured\n",
    "\n",
    "* We can choose the complexity of the model by changing the order of the polynomial model fit\n",
    "\n",
    "* Higher order polynomial models tend to have lower MSE on the data they're fit with\n",
    "\n",
    "---\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {}
   },
   "source": [
    "# Notation\n",
    "\\begin{align}\n",
    "x_{i} &\\quad \\text{input, independent variable}\\\\\n",
    "y_{i} &\\quad \\text{measurement, dependent variable}\\\\\n",
    "\\mathbf{x} &\\quad \\text{vector of input values}\\\\\n",
    "\\mathbf{y} &\\quad \\text{vector of measurements}\\\\\n",
    "\\hat{y}_{i} &\\quad \\text{estimate of dependent variable}\\\\\n",
    "\\epsilon_{i} &\\quad \\text{measurement error}\\\\\n",
    "\\theta &\\quad \\text{slope parameter}\\\\\n",
    "\\hat{\\theta} &\\quad \\text{estimated slope parameter}\\\\\n",
    "\\hat{\\theta}_\\text{MSE} &\\quad \\text{slope parameter estimated via the mean squared error}\\\\\n",
    "\\textrm{MSE} &\\quad \\text{mean squared error}\\\\\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "execution": {}
   },
   "source": [
    "---\n",
    "# Appendix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "execution": {}
   },
   "source": [
    "## Least Squares Optimization Derivation\n",
    "\n",
    "We will outline here the derivation of the least squares solution.\n",
    "\n",
    "We first set the derivative of the error expression with respect to $\\theta$ equal to zero, \n",
    "\n",
    "\\begin{align}\n",
    "\\frac{d}{d\\theta}\\frac{1}{N}\\sum_{i=1}^N(y_i - \\theta x_i)^2 = 0 \\\\\n",
    "\\frac{1}{N}\\sum_{i=1}^N-2x_i(y_i - \\theta x_i) = 0\n",
    "\\end{align}\n",
    "\n",
    "where we used the chain rule. Now solving for $\\theta$, we obtain an optimal value of:\n",
    "\n",
    "\\begin{align}\n",
    "\\hat\\theta = \\frac{\\sum_{i=1}^N x_i y_i}{\\sum_{i=1}^N x_i^2}\n",
    "\\end{align}\n",
    "\n",
    "Which we can write in vector notation as:\n",
    "\n",
    "\\begin{align}\n",
    "\\hat\\theta = \\frac{\\vec{x}^\\top \\vec{y}}{\\vec{x}^\\top \\vec{x}}\n",
    "\\end{align}\n",
    "\n",
    "\n",
    "This is known as solving the *normal equations*. For different ways of obtaining the solution, see the notes on [Least Squares Optimization](https://www.cns.nyu.edu/~eero/NOTES/leastSquares.pdf) by Eero Simoncelli."
   ]
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [],
   "include_colab_link": true,
   "name": "W4_Tutorial1",
   "provenance": [],
   "toc_visible": true
  },
  "interpreter": {
   "hash": "9516f62da91337f10c2adbe814d9c63a4b08f8271333386358218606edb781e3"
  },
  "kernel": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "kernelspec": {
   "display_name": "Python 3.7.11 64-bit ('pw3': conda)",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
